Synchronization Transition of Identical Phase Oscillators in a Directed Small- World 

Network 
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We numerically study a directed small-world network consisting of attractively coupled, identical 
phase oscillators. While complete synchronization is always stable, it is not always reachable from 
random initial conditions. Depending on the shortcut density and on the asymmetry of the phase 
coupling function, there exists a regime of persistent chaotic dynamics. By increasing the density of 
shortcuts or decreasing the asymmetry of the phase coupling function, we observe a discontinuous 
transition in the ability of the system to synchronize. Using a control technique, we identify the 
bifurcation scenario of the order parameter. We also discuss the relation between dynamics and 
topology and remark on the similarity of the synchronization transition to directed percolation. 



The adjustment of phase and frequency in large 
systems of oscillatory units can lead to global co- 
herent oscillations, i.e. synchronization. On the 
other hand, noise and heterogeneity in the system 
can weaken synchronization, or even destroy it. 
Synchronization in the nervous system can facili- 
tate the transfer of information or cause epileptic 
seizures. Multistability and hysteresis of normal 
and pathological collective behavior is observed. 
When all oscillators are identical and the cou- 
pling tends to decrease phase differences a state 
of complete synchronization is asymptotically sta- 
ble. But even in random networks with uncor- 
related and homogeneously distributed node de- 
grees this absorbing state may not be reached or 
disappear when it is perturbed locally. Here we 
perform a detailed numerical analysis of the tran- 
sition between different states of synchronization 
in a directed small-world network of phase os- 
cillators. By varying the mean in-degree of the 
network or the nonlinearity of the phase coupling 
function at zero phase difference, we find discon- 
tinuous and continuous transitions with mean- 
field critical behavior. 



I. Introduction 

Synchronization in spatially extended systems and 
complex networks is an important mechanism to create 
global spatio-temporal correlations from local interac- 
tion rules Its applications range from information 
transfer Q, self-organized optimization of work flow 
or traffic flow Q to the realization of strong coherent 
oscillations in arrays of Josephson junctions or coupled 
fiber lasers Q . The transition to partial synchronization, 
i.e., the emergence of a nonzero mean- field, in systems 
of nonidentical oscillators is well known. It has been 
studied analytically in the original texts by Kuramoto 
[1, Q and in a more general way in recent papers (9l-[ll| . 



Recently it has been shown that degree heterogeneity in 
scale-free random networks of identical oscillators can 
also lead to desynchronization or partially synchronized 
states even though complete synchronization is asymp- 
totically stable |12| . 

Here we study the transitions from incoherence to 
partial synchronization and to complete synchronization 
in a sparse, directed small-world network of identical 
phase oscillators. Since the formulation of the model jl^ 
small- world networks have been studied as a medium for 
dynamical processes, such as the spreading of epidemics 
p4| , the Ising model [Tsl - fTTj and also synchronization of 
nonidentical phase oscillators [l^. Nontrivial behavior 
in sparse networks of coupled dynamical systems is 
known to occur in boolean networks [l^, [1^ and has 
also been reported for leaky integrate-and-fire models 
with and without delay coupling [2l|, We show, 

for the case of attractively coupled, identical phase 
oscillators, that even if the in-degree distribution is 
homogeneous the system may enter a quasi-stationary 
chaotic state of incoherence or partial synchronization. 
So far, desynchronization, due to inherent heterogeneity, 
could be analytically described by global and local mean 
field approaches 0, which assume a sufficiently 

large number of neighbors. In contrast, our model is 
highly homogeneous, in terms of degree distribution and 
correlations, and none of the previously described mech- 
anisms for desynchronization in networks of oscillators 
can explain the observed transition. 



II. The Model 



The phase dynamics of N identical, weakly coupled 
limit-cycle oscillators may be described by the Kuramoto 
phase equations 0, 0, |1] . In a co- rotating frame of refer- 
ence where the common oscillator frequency is zero, these 
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FIG. 1: The network model. Unidirectional small-world net- 
works with A'^ = 32 nodes at (a), (b) low shortcut density 
(7 = 0.25 and (c), (d) higher shortcut density a — 2.0 . Joints 
of the network, i.e. nodes that receive more than one input, 
are marked gray, (b) At low shortcut densities most joints 
couple indirectly to two other joints through linear chain seg- 
ments of length L « 1/a. (d) At high shortcut densities each 
joint couples to fc = cr + 1 neighbors which are with high 
probability also joints. 



are 

N 

z?. = 5] i/., 5(^, -^0, (1) 

where Hij > is the coupling matrix and g{dj—'di) is the 
coupling function. The coupling function only depends 
on the phase difference between the oscillators and its 
shape is characteristic for given coupled limit-cycle os- 
cillators [3]. A typical coupling function, which is, for 
example, realized in diffusive coupling, has the proper- 
ties g(0) = and g'(0) > 0, implying that the coupling 
force between two coupled oscillators is attracting and 
vanishes when they have identical phases. For weakly an- 
harmonic oscillators the coupling function is usually well 
approximated by the first harmonics (See, e.g., (23|). i.e., 

5(At9) =sin(A?9-a)+sin(a). (2) 

The parameter a {—it/2 < a < 7r/2) breaks the anti- 
symmetry of the coupling function around zero, which 
can lead to nontrivial effects ^, [H, [23 - [26j . We can re- 
strict ourselves to nonnegative values of a (0 < a < 7r/2) 
since a transformation a — > —a and di — —di for all i 
leaves the Eqs. ([Ij and ^ invariant. 

Under very general conditions, i.e. nonnegative 
coupling strengths Hij > and a nondegenerate zero 
eigenvalue of the coupling Laplacian, complete synchro- 
nization with di — dj for all oscillators i and j is an 
asymptotically stable solution of Eq. ^ [23|. Strong 
connectedness of the network is a sufficient condition for 
this. Failure of a system of identical phase oscillators to 



synchronize can not be deduced from a local stability 
analysis of the completely synchronized state. It is 
the result of a chaotic phase dynamics that riddles the 
phase space so that complete synchronization cannot be 
reached. 

As a coupling network, we employ a well established 
modification JJ, 21, 2Sl of the original Watts Strogatz 
model [l^. Starting with a ring of N unidirectionally 
coupled oscillators we add A^sc unidirectional shortcuts 
with random origin i and destination j (See Fig. [!}. 
We refer to nodes which receive more than one input 
as the joints of the network. In addition to the system 
parameter a, which characterizes the phase coupling 
function, the shortcut density a = Ngc/N characterizes 
the topology of the coupling network. 

In the present paper, except for Fig. [2t>, we assume 
uniform and normalized input strength, i.e., 

N 

Y,H,, = l (3) 

for each oscillator and Hij = if oscillator i with 

in-degree fcj" couples to j or Hij = otherwise. The 
normalization would have a great impact on the dynam- 
ics if the coupling network was very heterogeneous. In 
the nonnormalized case and when the phases are uni- 
formly distributed the phase velocity of each oscillator is 
biased proportional to its in-degree and sin a [l^] ■ Oscil- 
lators with larger in-degree tend to move faster so that 
degree heterogeneity indirectly leads to a heterogeneity 
in frequencies. In contrast, when the coupling strength 
is normalized to unity the bias to the phase velocity is 
equal to sin a for all oscillators (See Fig. (^t). However, 
the directed small-world network model that we use has 
a Poissonian degree distribution, and is thus a homo- 
geneous network. In fact, we observe similar synchro- 
nization behavior with and without the normalization 
given by Eq. (See Figs. [2^ and[2}D). Another conse- 
quence of this normalization is that the phase velocities 
are bounded, i.e., —l<^i — sin a < 1 and an oscillator 
that receives only a single input is always phase locked 
to it. Phase slips do not occur in a chain of unidirection- 
ally coupled phase oscillators but only at the joints of a 
network. 



III. Simulation Results 

Figure [5] shows the phase diagrams obtained from 
numerical integrations of Eq. ([1} for networks of 
N = 800 oscillators starting from uniformly random 
initial conditions. For each set of parameter values for a 
and a, we show time and ensemble averages of various 
quantities. Numerical simulations of an ensemble of 
10 network realizations were run for T = 800 where 
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FIG. 2: Synchronization transition in the {a, a) parameter 
plane. Each point corresponds to an ensemble average over 
10 network realizations (A'' = 800) and time average over 600 
units after an initial transient of 200. Shown are (a) the mean 
order parameter R, (c) the mean oscillator frequency and (d) 
the variance of phase velocities in the case of normalized input 
strength. For comparison we also show (b) the mean order 
parameter for nonnormalized coupling strength for which a 
larger area of partial synchronization is observed at interme- 
diate shortcut densities. 



(a) 



a, 



. var(iJ)= 10" 



FIG. 3: Variance and distribution of phase velocities in the 
incoherent state, (a) Variance of phase velocities obtained 
from simulations with A'^ + Nsc = 10^ and a — -k jl (crosses) 
and from the simulations presented in Fig. [2ji (circles) at a = 
1.2. (b) The distribution of phase velocities in the incoherent 
state at (T = 1.0 (dots) is centered around the mean of = 
sina. It is peaked at the center and much broader than a 
Gaussian distribution of the same variance (dashed line) . 



an initial transient time of 200 time units was disre- 
garded. In Figs. [5^ and[5]D, the average order parameter 
R = (r-(t))timo,triais, where r{t) = | ^ e''''^*^ | is 
displayed. Note that r(t) = 1 for complete synchro- 
nization and i? '-^ 0{l/\fN) for complete incoherence, 
i.e., uniform phase distribution. There exists a clear 
transition from an incoherent regime {R ~ 0) at small 
shortcut densities a or larger asymmetry a to a coherent 
regime (i? = 1) at higher shortcut density or lower 
asymmetry. To quantify the dynamical properties of the 
incoherent state, we also observed the mean frequency 
{'&i{t))i, time, trials (Fig.Ht) and the variance of the phase 
velocities, (var ■di{t))i time, trials (Fig. [2jl). 

Figure ISb shows the distribution of phase veloci- 
ties, which is characteristic for the phase diffusion 




FIG. 4: Snapshots of the first 200 phases in a system of = 
1000 oscillators in a dynamical (quasi) equilibrium state for 
low shortcut density a = 0.35 (a — 0.36) in (a) and (b) 
and high shortcut density a = 128 (a = 1.5) in (c) and (d). 
Subfigures (a) and (c) show the stable incoherent state in a 
parameter region of bistability with the partially synchronized 
states shown in (b) and (d). Only at low shortcut densities 
the phases have a spatio-temporal structure at the length scale 
1/(7 of the chain segments. 



process in the incoherent state. The variance of the 
phase velocities is a measure for the internal noise due 
to chaotic phase dynamics. We observe that the mean 
frequency in the incoherent state with normalized input 
strength is equal to sina independently of a, while the 
variance of the phase velocities is independent of a. 

Figure S] shows snapshots of the first 200 phases in 
a network of = 1000 oscillators. In the incoherent 
state for low a (Fig. 0^), the one dimensional chain 
segments sustain traveling waves with an average phase 
difference of a between neighboring oscillators. Forward 
and backward phase slips occur occasionally in the joints 
of the network when the phase of the local mean field 
changes rapidly as it passes the vicinity of zero. In the 
partially synchronized state for low a (Fig. |4)d), the 
phases of the chain segments evolve around a global 
mean field. When a forward slip between the phase of a 
joint and the global mean field occurs, all oscillators in 
the chain segment adjacent to this joint will also perform 
a phase slip at a delayed time. For a high shortcut 
density no spatial patterns exist (Fig. |3J;, d). 



Control scheme and bifurcation scenario 

A more detailed examination using a linear control 
scheme reveals the bifurcation scenario for the order pa- 
rameter. Assuming that the order parameter r{t) evolves 
according to some unknown mean field dynamics, we 



FIG. 5: Bifurcation diagram of the order parameter i? as a 
function of control parameter q at selected values of shortcut 
densities a. Points on the branches of unstable (open squares) 
and stable (crosses) partially synchronized states were ob- 
tained as averages of the trajectory [R, a) (see light gray 
area in (d)) under the control scheme given by Eq. ((6|. The 
green lines are sixth order polynomial fits Q-{R) constrained 
to Q:'(0) = because of the assumption of a Hopf bifurca- 
tion of the incoherent state. From these fits we also find the 
threshold ap = a{l) for complete synchronization and the 
points of saddle node bifurcations asN of stable and unstable 
partially synchronized states. The dots in (a) are the average 
order parameter in simulation with networks of A'" = 8000 
oscillators. Each point is an ensemble average of the order 
parameter over 50 realizations after 1000 units of time. 



make the ansatz 



r = f{r,a,a,r]{t)), 



(4) 



where ri{t) represents intrinsic or external noise. In the 
absence of noise the dynamics under the general linear 
control scheme 



a 



co{r - To) + cif 



(5) 



has the fixed point r ~ ro and a — ao {ro , cr) with r = 
/(ro, q;o, cr, 0) = 0. In Appendix B we show that, in the 
absence of noise, sufficiently large positive values of cq 
and ci can stabilize any fixed point ro. In the presence 
of noise, the time derivative of r may not be well defined. 
We thus use a short-time average of r and f instead of 
their instantaneous values, i.e.. 



Co 



ro 



ci [r 



it)) 



7 ' 



(6) 



where 



(''(^))', = 7 / r{t — T)e ^'^dr, and 



(7) 



= 1 / f{t-T)e-''^dT 



FIG. 6: Numerically determined synchronization points for 

(a) low shortcut densities, (b) and (c) intermediate shortcut 
densities and (d) large shortcut densities. The open circles, 
the upward and the downward triangles mark the Hopf bi- 
furcation points q:hp(o") of the incoherent state, the transi- 
tion points ap (ct) to complete synchronization and the saddle- 
node bifurcation points asN(o"), respectively. Stable partial 
synchronization is found between ap{a) and asN(o"). At in- 
termediate to large shortcut densities (c) and (d) the transi- 
tion to synchronization is very well described by aHp(o") = 
arcsin((j/((T -I- 1)) (solid line). This is not the case for very 
low shortcut densities (a) where aup (cr) approaches zero more 
slowly than linearly. The line of slope 0.5 in the double- 
logarithmic plot (a) is drawn for comparison. The critical 
line q:hp(o") = arccos(0.85/Vo" + 1) obtained from the heuris- 
tic mean field ansatz Eq. (|13|l (dashed line in Subfig. (d)) 
agrees qualitatively with the asymptotic approach of q:hp to 
n/2 but is larger than the values obtained by our control 
scheme (open circles). The color code for the background of 

(b) and (c) is the same as in Fig. [2^ and b. 



are interpreted as stochastic integrals. 

Using the control scheme given by Eq. ([6]), and appro- 
priately tuned parameters 7, cq and ci, we succeeded to 
trace the stable and unstable branches in the bifurcation 
diagram of r as a function of the control parameter a 
at fixed shortcut densities a (Fig. [5]) . We find that the 
incoherent state always loses stability in a discontinuous 
transition at a critical value aHp(o'). We conjecture that 
this transition is a sub-critical Hopf-Bifurcation of the 
complex mean field. The branch of unstable partially 
synchronized states may fold back and become stable in 
a saddle-node bifurcation at a parameter asNlc). This 
hysteresis behavior is most pronounced around cr ss 1 
with possibly small parameter regions in which multiple 
partially synchronized states are stable (Fig. [St). At 
ap the branches of partial synchronization connect to 
the absorbing state of complete synchronization with 
r — 1. In the next section we perform a finite size scaling 
analysis at this transition point and conclude that it is 
of mean field directed percolation universality. From 
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these simulations, we could determine the three curves 
ckHp(o'), asN(o') and ap{a) for a large range of shortcut 
densities (Fig. [6]). 



Nonequilibrium Transition to Complete 
Synchronization 
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FIG. 7: Finite size scaling analysis of the nonequilibrium 
transition from partial to complete synchronization at ap 
for a — 1.0. In (a) simulation runs with system sizes up 
to A*' = 10*' were performed. At the critical point com- 
plete synchronization merges with the metastable partially 
synchronized state which is approached as if the tran- 
sition is of mean field directed percolation universality. The 
line (1 — r-(t)) ~ is drawn for comparison. The inset shows 
the linear approach of the mean order parameter r — >■ 1 in the 
vicinity of the critical point. The line of slope 1.01 is a lin- 
ear fit to the data in double- logarithmic scales. Other critical 
exponents are obtained from the time statistics for a realiza- 
tion to reach the absorbing state of complete synchronization. 
Figure (b) shows the median of this time at the critical point 
for various system sizes. In (c) we plot the fraction fj.o of 100 
realizations which reach complete synchronization before the 
time TN = 1900(Af/64000)° '''' for different as a function 
of a. This defines a median ai/2, shown in (d), which ap- 
proaches the critical point ap = 0.4065 at a power law with 
exponent —0.55 as a function of A'^. 

In finite systems, the partially synchronized state 
may disappear after a long transient time, and the 
state of complete synchronization is approached at 
an exponential rate. This sudden change of be- 
havior resembles the transition from an active state to 
the absorbing state in directed percolation processes [1^ . 

Finite size scaling analysis in low-dimensional cou- 
pling topologies has demonstrated that synchronization 
in coupled map lattices is of Kardar-Parisi-Zhang or 
directed percolation universality (30| . In a small world 
network, the transition is expected to be of mean field 
universality [isl [l8|. 



To examine this point we have performed a finite 
size scaling analysis in the vicinity of ap at fixed 
shortcut density a = 1.0 . The initial condition was 
a stable partially synchronized state at a = 0.52 and 
R ~ 0.75. After a transient time we decreased a and 
observed how r{t) approached zero. From simulation 
runs with single realizations of networks with sizes up 
to = 10^ oscillators in the vicinity of the transition 
point, we found ap — 0.4065 and (1 — i?) ~ (a — ap)'^ 
with /3 = 1.01 (Fig. [7]). This exponent is consistent 
with mean field universality of directed percolation. 
The nonequilibrium nature of the phase transition is 
expressed in the time dependence of the probability 
distribution of the order parameter. In a finite directed 
percolation system with a unique absorbing state, this 
absorbing state will be reached with probability one 
after a long enough transient time, even when a stable 
nonzero mean field solution exists 13 IP. 



A suitable absorbing state for a network of coupled 
phase oscillators can be defined via a Lyapunov function. 
Let : [0, 2tt]^ — [0, 2tt] denote the minimal length of 
the arc that contains the phases of all oscillators. Then 
one can show that is a Lyapunov function in any 
subset Bs of the phase space with V{Bs) < tt ~ 2a. 
If the network is strongly connected, i.e., there exists 
a directed path between any two nodes, and Hij > 
then one can show that y < for all ^ > in ^5. 
Thus Bs defines an absorbing state in which complete 
synchronization, i.e., = 0, is approached exponentially. 

Let fj.{a,N,t) denote the probability of a process 
to have reached the absorbing state prior to the time t, 
then the finite size scaling ansatz for this probability is 



/i(a, N,t)^fL ((a - ap)N'',N-H) 



(8) 



This probability can be determined by averaging over a 
large number of simulation runs. At each value of a in 
Fig. [T]: we have performed 100 simulations for system 
sizes 1000 < N < 64000. To determine the exponent z 
we have performed 1000 simulation runs for each system 
size at the critical parameter a = ap = 0.4065 and cr = 
1.0. At the transition point the median T1/2 defined as 
n{a,N,Ti/2) = 0.5 scales with the system size as 



(9) 



The third critical exponent v is related to the width and 
the position of the transition. Defining the position of 
the transition as ^{aif2, N^Tq) — 0.5 one finds 



(ai/2 - ap) N 



(10) 



Confidence intervals for the experimentally determined 
values of z and v could be estimated by bootstrapping on 
the sample of simulation runs. However, we presume that 
systematic errors are dominating for the relatively small 
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system sizes used in our simulations leading to slightly 
higher values of z and v than the expected value of 0.5 
for mean field directed percolation. We find z 0.54 and 
0.55. 



V. Discussion 

So far, we have presented our numerical results. In 
this section, we discuss the different dynamical regimes 
of our model and the mechanism for desynchronization. 




Topological crossover 

Our small-world network model exhibits a topological 
cross-over between low short cut densities a < 1 and 
high shortcut densities a > 1. The two regimes are 
characterized by the scaling of the average distance 
between successive joints of the network, which are oscil- 
lators that couple to more than one neighbor, along an 
arbitrary path (See Fig. [1]). By choosing the end points 
of the shortcuts randomly and uncorrelated, the number 
L of nodes between two joints on the original ring 
lattice is exponentially distributed as p{L) ~ exp(— ctL) 
in the limit iV — ^ oo. The expected value of L is 
(L) = (exp(cr) — 1) ^, which is approximately cr~^ at 
low shortcut densities and exp(— cr) at high shortcut 
densities. This cross-over is reflected in the phase 
dynamics of the incoherent state as a regime of traveling 
waves on the original ring lattice at cr < 1 to a regime 
without clear spatio-temporal patterns at ct > 1 (Fig. [4]). 

The Poissonian degree distribution with a mean in- 
degree of fc = cr + 1 and standard deviation ^/k provides 
homogeneity at high shortcut densities. On the other 
hand, at low shortcut densities cr < 1, the network con- 
sists of one dimensional chain segments which interact 
nonlinearly at the joints of the network. The network of 
joints, where chain segments are replaced by edges and 
indirect interactions, is also very homogeneous, since 
the number of joints which receive exactly two inputs is 
of order 0(cr) while the total number of all other joints 
is of order O(cr^). In this sense, the network is most 
inhomogeneous at intermediate shortcut densities cr sa 1 
where the amount of short chain segments and joints 
with more than two input links are comparable. Indeed, 
this topological complexity is reflected in a maximum 
of the variance of the phase velocities at intermediate 
shortcut densities, which quantifies the strength of 
the internal noise (See Fig. We will discuss the 

dynamical properties of our model for low and high 
shortcut densities, separately. 



FIG. 8: Complex correlation function in the incoherent state, 
(a) Absolute value in logarithmic scales as a function of the 
distance L on the ring backbone of the network, for k = a + 
1 = 4 (blue crosses), k — 6 (red circles) and k — 8 (green 
squares), (b) Angle as a function of the distance L on the 
ring backbone of the network for fc = 4 (blue crosses), fc = 6 
(red circles) and fc = 8 (green squares). We used a — 1.5, 
well above the synchronization threshold. The dashed lines 
in subfigures (a) and (b) mark the mean distance log N/ log fc 
in the network of size A'' = 32000. (c) Scaling of absolute 
value with the number of neighbors at distance L = 1 (blue 
crosses), L = 2 (red triangles) and L = 3 (green diamonds), 
(d) Logarithm of autocorrelation function at time difference 
r for fc = 25 (blue crosses), fc = 64 (red circles) and fc — 
100 (green triangles) and parametric fit to autocorrelation 
function of Brownian flight on the circle (Eq. (|ll|l . dashed 
lines). See Table I for values of D and 7. The inset shows the 
collapse of the curves under a rescaling of time with factor 
fc-1/2. 



Statistical properties of the incoherent state 

In the incoherent state, the phases undergo a chaotic 
phase diffusion process. An oscillator coupled to a finite 
number of neighbors is therefore subject to a fiuctuating 
force. In the incoherent state in a unidirectional network, 
the neighbors of an oscillator have independent phases, 
because their distance is of the order of the network di- 
ameter. Thus, if the number of neighbors k is sufficiently 
large, the local mean fields will fiuctuate around zero 
with approximately Gaussian distribution and one can 
show that its variance is equal to k~^. Since the ampli- 
tudes of the local mean fields determine the phase veloc- 
ities of the oscillators which in turn generate the local 
mean fields, the chaotic phase diffusion process is invari- 
ant under a rescaling of time as fc~^/^ (See Appendix A 
for details). From the circular autocorrelation function 
c(t) = (cos(i?(i + t) — d{t))) we can estimate the effec- 
tive phase diffusion constant D and an effective scattering 
rate 7 of the chaotic phase diffusion process by comparing 
it directly to the autocorrelation function of a Brownian 
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k 


25 


64 


100 




0.1252 


0.1296 


0.1417 




0.4064 


0.3937 


0.3598 


fcvar(i9) 


0.0478 


0.0519 


0.0544 


kDj = kv' 


0.0509 


0.0510 


0.0510 


CtHP 


1.27 


1.39 


1.43 



TABLE I: Timescales of the chaotic phase diffusion process for 
various large mean degrees fc ^ 1. We find that the effective 
phase diffusion constant and the effective scattering rate scale 
with l/\/fc and the variance of the phase velocities scales as 
1/k (See Fig. |3Jd). The transition point to synchronization 
ohp has been determined with the help of our control scheme. 



flight on the circle [32] with 
c(r) 



-_DT+D7"^(l-e""'^) 



(11) 



The variance of the phase velocities of such a Brownian 
flight is equal to = D7 and can directly be compared 
to the variance var('i$') of the chaotic phase diffusion 
process. In Table U we show the experimentally deter- 
mined effective diffusion constant and the variance of 
the phase velocities for k = 25, 64 and 100, and a = 7r/2. 
We find the predicted scaling of D, 7 and var(iS'). In 
particular the rescaled effective phase diffusion constant 
is DVk « 0.13. 

Note that these statistical properties do not depend on 
a as can be seen in Fig. 



Synchronization transition for high shortcut density 

In the incoherent state each oscillator in the net- 
work is subject to a fluctuating local mean field gen- 
erated by a finite number of neighbors performing the 
chaotic phase diffusion process. These neighbors sam- 
ple from the global distribution of phases. While the 
global mean field evolves deterministically in the thermo- 
dynamic limit N ^ 00, the local mean fields will be ap- 
proximately distributed as a Gaussian around the global 
mean field with a variance of order 0(l/fc), a relaxation 
rate 7 and a diffusion constant D of order 0{1/Vk). Let 
us define the complex correlation function 



Cij — {z^ Zj) 



(12) 



where Zj — exp(ii?j) and the average is over time. In our 
homogeneous network model, the correlation is a function 
cl of the length L of the shortest directed path from j to 
i. A heuristic mean field ansatz is to use this correlation 
function as a distance dependent weight and phase shift 
for the coupling 



R 



L = l 



(13) 



where R and O are the mean field amplitude and angle. 
Here we have chosen a co-rotating reference system where 
(i?) — 0. In the incoherent state the complex correlation 
function (Fig. [5]) has the form 



CL 



ALa 



(14) 



Each shell, i.e., a set of oscillators j with constant dis- 
tance L from i, is correlated to the oscillator i with an 
amplitude that decreases exponentially with L and at an 
angle which changes linearly with L. The factor go was 
determined to be approximately 0.85 (See Fig. [5t). We 
are now able to express the phase evolution as that of a 
globally coupled system with effective coupling strength 
Kcfi and effective asymmetry acs > a as 



= Koff R sin (9 - - tteff) , 
where k^s and are determined by 

(1 - r_-ia\L (1 - g)e' 



(15) 



E 

L = l 



ge 



1 ~ ge 



(16) 

This mean field ansatz predicts a stable incoherent state 
whenever > 7r/2 which is fullfilled if 



cos a < gok 2 . 



(17) 



The critical line anp — arccos^go / Vk) obtained from 
this heuristic mean field ansatz shows the same qual- 
itative behavior as the critical line obtained from the 
control scheme, although the numerical values do not 
agree well (See Fig. |6ji)- It predicts an asymptotic 
approach of a to 7r/2 at the order of 0(1/ Vk) as fc — ^ 00. 

Another heuristic argument can be made for the 
transition line ap{a). Suppose that all oscillators have 
the same phase = except one oscillator with phase 
which is externally driven. This oscillator assumes the 
role of an active seed in terms of percolation processes. 
If this active seed activates (i.e., induces oscillations of) 
the other oscillators that couple to the seed, the whole 
system may eventually become active, resulting in the 
incoherent or partially synchronized states. We are thus 
interested in the conditions under which an oscillator 
with phase and in-degree k that couples to the seed 
can become active. The phase equation for "di is initially 



1?! = ■isin(z9o 



^1 



k-l 



sin (z9i -I- a) -f 



sm a . 



(18) 

At criticality, we expect that the phase difference ??o — "^i 
is drifting fast so that the first term of Eq. (|18l) can be 
time-averaged and neglected. In this case, oscillator 1 
gets activated for sin a > (fc — 1)/A:, and with k — a + 1 
we obtain the condition 



a 



arcsm ■ 



(19) 
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for criticality. Equation ^T9\\ agrees surprisingly well with 
the numerically determined transition line q;hp(o') for 
a > 1 but less well with ap{a) (Fig. [Btjd). Again we 
predict that the critical a approaches 7r/2 at the order of 
0{l/Vk). 

Synchronization transition for low shortcut density 

For low shortcut densities cr < 1, the situation becomes 
even more complicated. Due to the topological cross-over 
another length scale is introduced in the system. 
The dynamics in the incoherent state for low shortcut 
densities is characterized by traveling waves along one 
dimensional chain segments of the original ring topology 
which interact nonlinearly at the joints of the network 
(See Fig. S]). We call the first and last oscillators in a 
chain segment the head and the tail of the chain, respec- 
tively. For normalized input strength the chain segments 
are always phase locked to the head but act as a low 
pass filter to the phase dynamics (Appendix C) . As a re- 
sult, the variance of the phase velocities decreases along 
a chain segment and for smaller values of a (Fig.[3l). The 
interaction between the joints of the network is indirect 
and involves delay and additional phase shifts. Also, the 
correlation between the head and an oscillator down the 
chain segment decreases exponentially with the distance 
between the head and the oscillator. These factors seem 
to inhibit synchronization so that the critical q;//p((t) is 
still decreasing for cr <C 1 but much slower than linearly 
(Fig.Eti) . 

VI. Summary and Conclusions 

We have shown that a finite average number of neigh- 
bors and an asymmetry of the phase coupling function 
can inhibit synchronization in homogeneous networks 
of identical oscillators. Using a control technique, we 



could construct the bifurcation diagram for the order 
parameter. 

A Finite size scaling analysis at the nonequilibrium 
phase transition from partial synchronization to com- 
plete synchronization shows critical exponents of mean 
field universality. 

In contrast to it is not the heterogeneity of 

the node degree distribution that drives the system 
away from synchronization. Instead, it is the inter- 
play between the spatial structure of the network and 
the internal noise which prevents the oscillators from 
reaching synchronization. The temporal fluctuations 
generated by the system itself in the chaotic incoherent 
or partially synchronized state give rise to a correlation 
function that decays exponentially with the distance. 
Using this correlation function to formulate a heuristic 
mean-field theory for the incoherent state, we have 
qualitatively explained the transition from incoherence 
to synchronization. 

Our main results are derived from numerical simu- 
lations. An analytic description of the transition curve 
and the chaotic phase diffusion process in the incoherent 
state remain challenging open problems. 
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Appendix A : Scaling Argument for ct ^ 1 

Here we will present the details of the scaling argument 
that relates the amplitude and the time-scale of the fluc- 
tuations in the local fields. Let us consider the dynamics 
of the complex state Zi — exp(i?9i) of an oscillator 



i Im |z*^ie 



r, i Im 



where '^i is the local field given by 



(20) 
(21) 



where (...) indicates the time average with respect to 
the stationary processes Zj{t). Assuming k independent 
neighbors with correlation function 



c,(r) = (z*(i + r)z,(t)) 



one finds exactly 



(23) 



(24) 



This relation is remarkable since it is valid for all k and 
all time differences r. Therefore, the phase dynamics 
of the local fields has exactly the same time scale 
as the dynamics of the phase oscillators. Without 
calling on the central limit theorem this relation also 
gives the mean square amplitude of the local fields as 
(r2) = c*(0) = k-^ for ah k. 

The assumption of independent phases of the neighbors 
only holds in large unidirectional networks in the inco- 
herent state. Due to the chaotic dynamics, correlations 
between oscillators decay exponentially with the distance 
in the network which is of order 0(log N/ log k) (Fig. [5]). 



The effective phase diffusion constant 
fined by the asymptotic behavior of c*(t) 
large r 



can be de- 
and Cz{t) at 



D 



lim — In c(r) 



(25) 



We recall that the circular autocorrelation function 
for a free Brownian flight on the circle is c(r) = 
exp {—Dt + D/j (1 — exp(— 7r))), where 7 is an effective 
scattering rate and D is the effective phase diffusion con- 
stant |32l] . Exponential decay is expected at time scales 
T7 » 1. From Eqs. and (US]) it follows that the ef- 
fective phase diffusion constants of the single oscillator 
phase and that of the local flelds are identical. For suf- 
ficiently large fc, the local fields are normally distributed 
and it is only the amplitude of the local fields that de- 
creases when k is increased. The process (|20|) is therefore 
invariant under a rescaling 



t' 



*'(t') 



k 

k^' 



(26) 



In our model, this rescaling is achieved by changing 
the number of neighbors from k to k' . The second 
transformation affects both the amplitude and the phase 
diffusion timescale of the local fields. 



The amplitude and the time scale of local field dynamics 
can be inferred from the complex correlation function 



C*(t) = (vI/*(t + T)vI/,(t)) , 



(22) 



Because of the homogeneity of our model, we can 
assume statistical equivalence of the phase dynamics for 
all oscillators. The phase 'di{t) of an oscillator that is 
coupled via Eq. (PH)) to k independent but identically 
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FIG. 9: Effective phase difTusion constant Dz of a pfiase 
oscillator with complex state variable z coupled to a com- 
plex valued Ornstein-Uhlenbeck process of unit variance and 
phase diffusion constant (Eq. (I27|l 'l. The fixed point at 
D* ~ 0.15 is expected to be close to the rescaled effective 
phase diffusion constant D{k)\/k of the characteristic station- 
ary phase diffusion process in the incoherent state for fc ^ 1. 



distributed phase diffusion processes must be 

another reahzation of the same process. In particular, 
the effective phase diffusion constant must be the same. 
To get insight into the phase diffusion process, let us 
examine an oscillator with complex state z that is 
coupled to a complex Ornstein-Uhlenbeck process 4' 



= izlm[z*\l'] 



(27) 



where rj is complex Gaussian white noise of unit 
strength in real and imaginary part. Note that the 
Ornstein-Uhlenbeck process ^ is only an approximation 
of the dynamics of an actual local field which has the 
same Gaussian stationary distribution with expected 
square amplitude (r'^) — 1. The effective phase diffusion 
constant as defined in Eqs. (|22p and (l^5|) is = D 
[3^ . The effective phase diffusion constant depends 
nonlinearly on . For large values of D the time scales 
of ^ and z are separated and the phase of z will diffuse 
much slower than the phase of For very small values 
of D the phase of z will almost always be locked to 
the phase of Only when 4' diffuses near zero, phase 
slips may occur which add a ballistic component to 
the evolution of z and can increase the effective phase 
diffusion constant to a value larger than D. 

There exists a fixed point D* for which the time 
scale of ^' and z are identical. The fixed point D* 
should be close to the actual value of the self-consistent 
solution. A self-consistent solution must also yield 



Eq. (|24p . The average field of k such processes z may be 
normally distributed but is not an Ornstein-Uhlenbeck 
process. However, the iterative procedure of coupling 
a phase oscillator to k independent neighbors defines 
a mapping in the space of stationary processes on the 
circle. Starting with a Brownian phase diffusion we see 
that in the first iteration the map is contracting with 
respect to the effective diffusion constant (Fig. [S]), an 
indication for the existence of a unique stationary phase 
diffusion process characterizing the incoherent state. 

We have measured as a function of D in nu- 
merical integration of Eq. ([27)) and found the the fixed 
point to be D* w 0.15. Table U shows the effective 
diffusion constants in our model for k = 25, 64 and 100, 
obtained from a nonlinear parametric least square fit to 
the autocorrelation function (system average). We find 
that Dy/k « 0.13 agrees well with this value. 



Appendix B : Control Scheme 



Given the mean field dynamics 
r = f{r,a) 



(28) 



of a scalar order parameter r with a system parameter 
a we are interested in the bifurcation curve ao{r) with 
/(r, ao(r)) = 0. If we can control the system parameter 
a = a{t) then the linear control scheme 



a 



co{r - To) + cif 



(29) 



has the fixed point (r, a) — (ro,ao(ro)). With /o — 
drfiro,Oio{ro)) and /i = daf{ro,ao{ro)) the Jaccobian 
of the combined mean field and control dynamics reads 



J 




(30) 



Assuming fi < 0, the conditions for stability tr(J) < 
and Det(J) > yield 



ci > and Co > /oci. 

Ji 



(31) 



The assumption of /i < applies, since an increase of a 
counteracts synchronization and decreases r if a > aQ^r). 
Sufficient conditions to stabilize any point (ro, ao(''o)), 
regardless of the sign of /o, i.e., of the stability of the 
uncontrolled fixed point, are 



CO > |/o|ci > f^. 



1 
l/ll 



(32) 



Appendix C : The overdamped linear chain 

For low shortcut densities one can view the dynamics 
in the small-world network as traveling waves which 
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interact through a network of joints. Each joint of the 
network is the head of a unidirectional chain segment. 
A joint receives input from at least two nodes of other 
chain segments. To understand the role of the chain 
segments in the transition to synchronization, we study 
the dynamics of the phases and the signal transmission 
along a chain segment in a linear approximation. 



The phase equations for 
chain of oscillators are 



unidirectionally coupled 



= sin (z?„_i - i?„ - a) + 



sma 



(33) 



We identify ??o with the phase in the head oscillator of the 
chain and with the phase of the tail oscillator. Un- 
der appropriate boundary conditions Eq. (|33)) allows for 
traveling wave solutions d = sin a and i9„_i — = a. In 
principle, all frequencies in the interval [sina— 1, sina+l] 
and corresponding phase differences are possible but we 
choose the traveling wave solution with the average fre- 
quency sin a. Small deviations Xn from this solution can 
be studied in a linear approximation 



^n— 1 



(34) 



which can be solved given the trajectory XQ{t) in the head 
of the chain. The system is asymptotically independent 
from initial conditions and the solution is 



Xn{t) ^ / e '^Xn-l{t~T)dT ^ / t" ^ — XQ{t~T)dT 



(n-iy. 



(35) 

The dynamics in the tail of a chain is a time convolution 
of the dynamics in the head of the same chain around 
the traveling wave solution. A discontinuous jump of the 
phase in the head of the chain will translate at unit speed 
while the width of the phase jump grows at the same 
rate, as a result of the gamma distribution for the time 
convolution kernel. In this linear approximation, the de- 
viation from the traveling wave solution grows diffusively 
with the distance from the head of the chain. Suppose 
the head of the chain performs a Brownian motion with 
{(xoit) - xo(0))2) = 2Dt. Then we obtain 



and 



= / dTdT'^ ' 





r2(n) 

X (:,o(-T).To(-T'))|,„(o)=o 



OO 

■J- 



2 J dTdr 





iT{T + T)y 



r2(r 



^ (^o( ^))|:eo(o)=0 

OO 

1 (n) 1 [n) 



Using 



r(n+ l,r) = nr(n,T) -fT"e 



n-T 



(37) 



(38) 



we obtain 



|a;o(0)=0 



AD J dT 



X ' 



r(n + i)r(n + i) 
(r(n + i,r) -r"e"^) 



2Dn 1-2 



2n„-2T 



r2(n-|- 1) 



2Dn 1 - 



r2(n + 1) 

1 



(39) 



The effective spatial phase diffusion constant along a lin- 
ear chain of oscillators is the same as the temporal diffu- 
sion constant of the phase in the head of the chain. The 
complex correlation function 



CLo(O) = {zUt)zoit)) 



JaL-DL 



for L > 1 (40) 



has an amplitude that decreases exponentially with the 
chain length. A lower shortcut density therefore de- 
creases the correlation between the joints of the network, 
making it more difficult to synchronize. 



{{Xnit) - Xoit)f) = (Xnitf) 



xo{t)=0 ■■ 



(36) 



